library(tidyverse)
library(lme4)
library(lmerTest)
library(knitr)
library(cowplot)
library(caret)
library(ROCR)algorithm = c("#006989", "#FEC601", "#F43C13", "#00A5CF", "#00A878")
instruction = wesanderson::wes_palette("Darjeeling1", 2, "continuous")
craving = wesanderson::wes_palette("Darjeeling1", 3, "continuous")
rating = c("#00A08A", "#F2AD00", "#F98400", "#FF0000")
dc_bw = plot_aes = theme_minimal() +
theme(legend.position = "top",
legend.text = element_text(size = 12),
text = element_text(size = 16, family = "Futura Medium"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text = element_text(color = "black"),
axis.line = element_line(colour = "black"),
axis.ticks.y = element_blank())source("load_data.R")Check if wrong buttons were used (i.e., not 5-8)
NB! NEED TO GO THROUGH NEW PARTICPANTS/WAVES
subs = data.all %>%
group_by(subjectID, wave, run, rating) %>%
summarize(n = n()) %>%
spread(rating, n) %>%
mutate(messed = ifelse(is.na(`5`) & !is.na(`<NA>`), "yes", NA)) %>%
filter(messed == "yes") %>%
ungroup() %>%
select(subjectID, wave) %>%
unique()
data.all %>%
group_by(subjectID, run, rating) %>%
summarize(n = n()) %>%
spread(rating, n) %>%
mutate(messed = ifelse(is.na(`5`) & !is.na(`<NA>`), "yes", NA)) %>%
filter(subjectID %in% subs$subjectID)Recoding
* DEV069: recode runs1
NB! NEED TO GO THROUGH NEW PARTICPANTS/WAVES
data.ex = data.all %>%
mutate(rating = ifelse(subjectID == "DEV069" & run == "run1", rating - 1, rating),
rating = ifelse(subjectID == "DEV069" & run == "run1" & is.na(rating), 8, rating),
rating = rating - 4) %>%
group_by(subjectID, wave) %>%
arrange(subjectID, run) %>%
mutate(trial = row_number())file_dir = "~/Documents/code/sanlab/DEV_scripts/fMRI/betaseries/ROC/dotProducts_ROC_wave1/"
file_pattern = "DEV[0-9]{3}_meanIntensity.txt"
file_list = list.files(file_dir, pattern = file_pattern)
intensities = data.frame()
for (file in file_list) {
temp = tryCatch(read.table(file.path(file_dir,file), fill = TRUE) %>%
rename("subjectID" = V1,
"meanIntensity" = V3) %>%
extract(V2, "beta", "beta_([0-9]{4}).nii") %>%
mutate(beta = as.integer(beta),
wave = 1), error = function(e) message(file))
intensities = rbind(intensities, temp)
rm(temp)
}
file_dir = "~/Documents/code/sanlab/DEV_scripts/fMRI/betaseries/ROC/dotProducts_ROC_wave2/"
file_list = list.files(file_dir, pattern = file_pattern)
for (file in file_list) {
temp = tryCatch(read.table(file.path(file_dir,file), fill = TRUE) %>%
rename("subjectID" = V1,
"meanIntensity" = V3) %>%
extract(V2, "beta", "beta_([0-9]{4}).nii") %>%
mutate(beta = as.integer(beta),
wave = 2), error = function(e) message(file))
intensities = rbind(intensities, temp)
rm(temp)
}file_dir = "~/Documents/code/sanlab/DEV_scripts/fMRI/betaseries/ROC/dotProducts_ROC_wave1/"
file_pattern = "DEV[0-9]{3}_dotProducts.txt"
file_list = list.files(file_dir, pattern = file_pattern)
dots = data.frame()
for (file in file_list) {
temp = tryCatch(read.table(file.path(file_dir,file), fill = TRUE) %>%
rename("subjectID" = V1,
"map" = V3,
"dotProduct" = V4) %>%
extract(V2, "beta", "beta_([0-9]{4}).nii") %>%
extract(map, "algorithm", "(.*)_.*.nii") %>%
mutate(beta = as.integer(beta),
wave = 1), error = function(e) message(file))
dots = rbind(dots, temp)
rm(temp)
}
file_dir = "~/Documents/code/sanlab/DEV_scripts/fMRI/betaseries/ROC/dotProducts_ROC_wave2/"
file_list = list.files(file_dir, pattern = file_pattern)
for (file in file_list) {
temp = tryCatch(read.table(file.path(file_dir,file), fill = TRUE) %>%
rename("subjectID" = V1,
"map" = V3,
"dotProduct" = V4) %>%
extract(V2, "beta", "beta_([0-9]{4}).nii") %>%
extract(map, "algorithm", "(.*)_.*.nii") %>%
mutate(beta = as.integer(beta),
wave = 2), error = function(e) message(file))
dots = rbind(dots, temp)
rm(temp)
}dots.merged = dots %>%
left_join(., intensities, by = c("subjectID", "wave", "beta")) %>%
group_by(subjectID, wave, algorithm) %>%
mutate(rownum = row_number())
# plot original
dots.merged %>%
filter(algorithm == "craving_regulation") %>%
ggplot(aes(1, meanIntensity)) +
geom_boxplot()# assess extreme values and exclude when calculating SDs
dots.merged %>%
filter(algorithm == "craving_regulation") %>%
arrange(meanIntensity)dots.merged %>%
filter(algorithm == "craving_regulation") %>%
arrange(-meanIntensity)# recode outliers as NA
dots.merged = dots.merged %>%
ungroup() %>%
mutate(meanIntensity = ifelse(meanIntensity > 1 | meanIntensity < -1, NA, meanIntensity),
median = median(meanIntensity, na.rm = TRUE),
sd3 = 3*sd(meanIntensity, na.rm = TRUE),
outlier = ifelse(meanIntensity > median + sd3 | meanIntensity < median - sd3, "yes", "no"),
dotProduct = ifelse(outlier == "yes", NA, dotProduct))
# plot after
dots.merged %>%
filter(algorithm == "craving_regulation") %>%
ggplot(aes(1, meanIntensity)) +
geom_boxplot()NB! NEED TO GO THROUGH NEW PARTICPANTS/WAVES
trial.numbers = data.frame(subjectID = c(rep("DEV060", 79), rep("DEV061", 79), rep("DEV063", 71), rep("DEV081", 80), rep("DEV082", 75)),
wave = 1,
rownum = c(1:79, 1:79, 1:71, 1:80, 1:75),
trial = c(1:19, 21:80, 1:59, 61:80, 1:31, 41:80, 1:20, 41:80, 21:40, 1:35, 41:80))
dots.check = dots.merged %>%
group_by(subjectID, wave, algorithm) %>%
mutate(rownum = row_number()) %>%
left_join(., trial.numbers, by = c("subjectID", "wave", "rownum")) %>%
mutate(trial = ifelse(is.na(trial), rownum, trial),
dotProduct = ifelse(subjectID == "DEV060" & wave == 1 & trial %in% 19:20, NA,
ifelse(subjectID == "DEV061" & wave == 1 & trial == 59, NA,
ifelse(subjectID == "DEV082" & wave == 1 & trial %in% 19:20, NA, dotProduct)))) %>%
select(-rownum) #%>%
#left_join(., striping, by = c("subjectID", "beta")) %>%
#mutate(dotProduct = ifelse(!is.na(striping), NA, dotProduct))dots.ex = dots.check %>%
group_by(subjectID, algorithm) %>% # standardize within sub and algorithm
mutate(dotSTD = scale(dotProduct, center = FALSE)) Exclusions
Other
* select only craved trials
NB! NEED TO GO THROUGH NEW PARTICPANTS/WAVES
data_all = left_join(dots.ex, data.ex, by = c("subjectID", "wave", "trial")) %>%
filter(!(subjectID %in% c("DEV001","DEV020","DEV032","DEV047","DEV055","DEV064","DEV066", "DEV017", "DEV019", "DEV037", "DEV054") & wave == 1)) %>%
filter(!(subjectID == "DEV029" & wave == 1 & run == "run3") & !(subjectID == "DEV037" & wave == 1 & run == "run1") & !(subjectID == "DEV042" & wave == 1 & run == "run4") & !(subjectID == "DEV067" & wave == 1 & run == "run4")) %>%
ungroup() %>%
mutate(algorithm = gsub("_signature", "", algorithm),
wave = as.character(wave))
data = data_all %>%
filter(craving == "craved")data %>%
select(subjectID, wave) %>%
unique() %>%
group_by(wave) %>%
summarize(n = n())data %>%
filter(algorithm == "craving_regulation") %>%
group_by(subjectID) %>%
summarize(n = n()) %>%
arrange(n)# roc curve
perf.df = data %>%
filter(algorithm == "craving_regulation") %>%
filter(!is.na(dotSTD) & !is.na(craving)) %>%
mutate(instruction = ifelse(instruction == "regulate", 1, 0)) %>%
group_by(wave) %>%
do({
wave = .$wave
pred = prediction(.$dotSTD, .$instruction, label.ordering = NULL)
perf = performance(pred, measure = "tpr", x.measure = "fpr")
data.frame(cut=perf@alpha.values[[1]],fpr=perf@x.values[[1]],tpr=perf@y.values[[1]])
})
ggplot(perf.df, aes(fpr, tpr, color = as.factor(wave))) +
geom_line() +
scale_color_manual(values = algorithm) +
xlab("false positive rate") +
ylab("true positive rate")roc = data %>%
filter(algorithm == "craving_regulation") %>%
select(subjectID, trial, instruction, rating, dotSTD) %>%
mutate(guess.instruction = ifelse(dotSTD > 0, "regulate", "look"),
guess.rating = ifelse(dotSTD > 0, "low", "high"),
rating.bin = ifelse(rating >= 3, "high", "low"),
instruction = as.factor(instruction),
guess.instruction = as.factor(guess.instruction),
rating.bin = as.factor(rating.bin),
guess.rating = as.factor(guess.rating))
confusionMatrix(roc$guess.instruction, roc$instruction)## Confusion Matrix and Statistics
##
## Reference
## Prediction look regulate
## look 5551 2292
## regulate 3490 6707
##
## Accuracy : 0.6795
## 95% CI : (0.6726, 0.6863)
## No Information Rate : 0.5012
## P-Value [Acc > NIR] : < 0.00000000000000022
##
## Kappa : 0.3592
##
## Mcnemar's Test P-Value : < 0.00000000000000022
##
## Sensitivity : 0.6140
## Specificity : 0.7453
## Pos Pred Value : 0.7078
## Neg Pred Value : 0.6577
## Prevalence : 0.5012
## Detection Rate : 0.3077
## Detection Prevalence : 0.4348
## Balanced Accuracy : 0.6796
##
## 'Positive' Class : look
##
#confusionMatrix(roc$guess.rating, roc$rating.bin)# roc curve
perf.df = data %>%
filter(algorithm == "craving") %>%
filter(!is.na(dotSTD) & !is.na(craving)) %>%
mutate(instruction = ifelse(instruction == "regulate", 1, 0)) %>%
group_by(wave) %>%
do({
wave = .$wave
pred = prediction(.$dotSTD, .$instruction, label.ordering = NULL)
perf = performance(pred, measure = "tpr", x.measure = "fpr")
data.frame(cut=perf@alpha.values[[1]],fpr=perf@x.values[[1]],tpr=perf@y.values[[1]])
})
ggplot(perf.df, aes(fpr, tpr, color = as.factor(wave))) +
geom_line() +
scale_color_manual(values = algorithm) +
xlab("false positive rate") +
ylab("true positive rate")roc = data %>%
filter(algorithm == "craving") %>%
select(subjectID, trial, instruction, rating, dotSTD) %>%
mutate(guess.instruction = ifelse(dotSTD > 0, "regulate", "look"),
guess.rating = ifelse(dotSTD > 0, "low", "high"),
rating.bin = ifelse(rating >= 3, "high", "low"),
instruction = as.factor(instruction),
guess.instruction = as.factor(guess.instruction),
rating.bin = as.factor(rating.bin),
guess.rating = as.factor(guess.rating))
confusionMatrix(roc$guess.instruction, roc$instruction)## Confusion Matrix and Statistics
##
## Reference
## Prediction look regulate
## look 3629 4178
## regulate 5412 4821
##
## Accuracy : 0.4684
## 95% CI : (0.4611, 0.4757)
## No Information Rate : 0.5012
## P-Value [Acc > NIR] : 1
##
## Kappa : -0.0629
##
## Mcnemar's Test P-Value : <0.0000000000000002
##
## Sensitivity : 0.4014
## Specificity : 0.5357
## Pos Pred Value : 0.4648
## Neg Pred Value : 0.4711
## Prevalence : 0.5012
## Detection Rate : 0.2012
## Detection Prevalence : 0.4328
## Balanced Accuracy : 0.4686
##
## 'Positive' Class : look
##
#confusionMatrix(roc$guess.rating, roc$rating.bin)# roc curve
perf.df = data_all %>%
filter(algorithm == "craving") %>%
filter(!is.na(dotSTD)) %>%
filter(craving %in% c("neutral", "craved")) %>%
mutate(craving = ifelse(craving == "craved", 1, 0)) %>%
group_by(wave) %>%
do({
wave = .$wave
pred = prediction(.$dotSTD, .$craving, label.ordering = NULL)
perf = performance(pred, measure = "tpr", x.measure = "fpr")
data.frame(cut=perf@alpha.values[[1]],fpr=perf@x.values[[1]],tpr=perf@y.values[[1]])
})
ggplot(perf.df, aes(fpr, tpr, color = as.factor(wave))) +
geom_line() +
scale_color_manual(values = algorithm) +
xlab("false positive rate") +
ylab("true positive rate")roc = data_all %>%
filter(algorithm == "craving") %>%
filter(craving %in% c("neutral", "craved")) %>%
select(subjectID, trial, craving, rating, dotSTD) %>%
mutate(guess.craving = ifelse(dotSTD > 0, "craved", "neutral"),
guess.rating = ifelse(dotSTD > 0, "low", "high"),
rating.bin = ifelse(rating >= 3, "high", "low"),
craving = as.factor(craving),
guess.craving = as.factor(guess.craving),
rating.bin = as.factor(rating.bin),
guess.rating = as.factor(guess.rating))
confusionMatrix(roc$guess.craving, roc$craving)## Confusion Matrix and Statistics
##
## Reference
## Prediction craved neutral
## craved 10233 3897
## neutral 7807 5143
##
## Accuracy : 0.5678
## 95% CI : (0.5619, 0.5737)
## No Information Rate : 0.6662
## P-Value [Acc > NIR] : 1
##
## Kappa : 0.1229
##
## Mcnemar's Test P-Value : <0.0000000000000002
##
## Sensitivity : 0.5672
## Specificity : 0.5689
## Pos Pred Value : 0.7242
## Neg Pred Value : 0.3971
## Prevalence : 0.6662
## Detection Rate : 0.3779
## Detection Prevalence : 0.5218
## Balanced Accuracy : 0.5681
##
## 'Positive' Class : craved
##
#confusionMatrix(roc$guess.rating, roc$rating.bin)data %>%
filter(!is.na(rating)) %>%
ggplot(aes(algorithm, dotSTD, fill = instruction)) +
stat_summary(fun.y = mean, geom = "bar", position = position_dodge(width = 0.95)) +
stat_summary(fun.data = mean_cl_boot, geom = "errorbar", position = position_dodge(width = 0.95), width = 0) +
scale_fill_manual(name = "", values = instruction) +
labs(y = "standardized pattern expression value\n", x = "") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
ggplot(aes(instruction, dotSTD)) +
stat_summary(aes(group = subjectID), fun.y = mean, geom = "line", alpha = .1, size = .5) +
stat_summary(aes(group = 1), fun.y = mean, geom = "line", size = .75) +
stat_summary(aes(color = instruction), fun.data = "mean_cl_boot", geom = "pointrange", width = 0, size = .75) +
facet_grid(~algorithm) +
scale_color_manual(name = "", values = instruction) +
labs(y = "standardized pattern expression value\n", x = "") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
ggplot(aes(instruction, dotSTD, fill = wave)) +
stat_summary(fun.y = mean, geom = "bar", position = position_dodge(width = 0.95)) +
stat_summary(fun.data = mean_cl_boot, geom = "errorbar", position = position_dodge(width = 0.95), width = 0) +
facet_grid(~algorithm) +
scale_fill_manual(name = "wave", values = algorithm) +
labs(y = "standardized pattern expression value\n", x = "") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
ggplot(aes(wave, dotSTD, color = instruction)) +
stat_summary(aes(group = interaction(subjectID, instruction)), fun.y = mean, geom = "line", alpha = .1, size = .5) +
stat_summary(aes(group = instruction), fun.y = mean, geom = "line", size = .75) +
stat_summary(fun.data = "mean_cl_boot", geom = "pointrange", width = 0, size = .75) +
facet_grid(~algorithm) +
scale_color_manual(name = "", values = instruction) +
labs(y = "standardized pattern expression value\n", x = "") +
dc_bwdata %>%
ggplot(aes(dotSTD, rating, color = algorithm, fill = algorithm)) +
geom_jitter(alpha = .1, height = .1) +
geom_smooth(method = "lm") +
scale_color_manual(name = "signature", values = algorithm) +
scale_fill_manual(name = "signature", values = algorithm) +
labs(y = "craving rating\n", x = "\nstandardized pattern expression value") +
dc_bwdata %>%
ggplot(aes(dotSTD, rating, color = wave, fill = wave)) +
geom_jitter(alpha = .1, height = .1) +
geom_smooth(method = "lm") +
facet_grid(~algorithm) +
scale_color_manual(name = "wave", values = algorithm) +
scale_fill_manual(name = "wave", values = algorithm) +
labs(y = "craving rating\n", x = "\nstandardized pattern expression value") +
dc_bwdata %>%
ggplot(aes(dotSTD, rating, color = instruction, fill = instruction)) +
geom_jitter(alpha = .1, height = .1) +
geom_smooth(method = "lm") +
facet_grid(~algorithm) +
scale_color_manual(name = "instruction", values = instruction) +
scale_fill_manual(name = "instruction", values = instruction) +
labs(y = "craving rating\n", x = "\nstandardized pattern expression value") +
dc_bwdata %>%
ggplot(aes(dotSTD, rating, color = wave, fill = wave)) +
geom_jitter(alpha = .05, height = .1) +
geom_smooth(method = "lm") +
facet_grid(instruction~algorithm) +
scale_color_manual(name = "wave", values = algorithm) +
scale_fill_manual(name = "wave", values = algorithm) +
labs(y = "craving rating\n", x = "\nstandardized pattern expression value") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
ggplot(aes(dotSTD, rt, color = instruction, fill = instruction)) +
geom_jitter(alpha = .05, height = .1) +
geom_smooth(method = "lm", alpha = .2) +
facet_grid(~algorithm) +
scale_color_manual(name = "", values = instruction) +
scale_fill_manual(name = "", values = instruction) +
labs(x = "\nstandardized pattern expression value", y = "reaction time (seconds)\n") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
filter(algorithm == "craving_regulation") %>%
ggplot(aes(rating, rt, color = instruction, fill = instruction)) +
geom_jitter(alpha = .05, height = .1) +
geom_smooth(method = "lm", alpha = .2) +
facet_grid(~algorithm) +
scale_color_manual(name = "", values = instruction) +
scale_fill_manual(name = "", values = instruction) +
labs(x = "\ncraving rating", y = "reaction time (seconds)\n") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
ggplot(aes(dotSTD, rt, color = wave, fill = wave)) +
geom_jitter(alpha = .05, height = .1) +
geom_smooth(method = "lm", alpha = .2) +
facet_grid(instruction~algorithm) +
scale_color_manual(name = "", values = algorithm) +
scale_fill_manual(name = "", values = algorithm) +
labs(x = "\nstandardized pattern expression value", y = "reaction time (seconds)\n") +
dc_bwdata %>%
filter(!is.na(rating)) %>%
filter(algorithm == "craving_regulation") %>%
ggplot(aes(rating, rt, color = wave, fill = wave)) +
geom_jitter(alpha = .05, height = .1) +
geom_smooth(method = "lm", alpha = .2) +
facet_grid(~instruction) +
scale_color_manual(name = "", values = algorithm) +
scale_fill_manual(name = "", values = algorithm) +
labs(x = "\ncraving rating", y = "reaction time (seconds)\n") +
dc_bwWhen looking, people showed increases in expression of the craving regulation signature (i.e., their brains looks more like they’re regulating) from wave 1 to 2.
When regulating, people showed decreases in expression of the craving regulation signature (i.e., their brains look less like they’re regulating) from wave 1 to 2. This may mean that it’s gotten easier and/or they’re exprience less craving and therefore need to regulate less.
mod_pev_wave = lmer(dotSTD ~ wave * instruction + (1 + wave * instruction | subjectID),
data = filter(data, algorithm == "craving_regulation"))
ggeffects::ggpredict(mod_pev_wave, terms = c("wave", "instruction")) %>%
data.frame() %>%
ggplot(aes(x, predicted, color = group, fill = group)) +
geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
geom_line(aes(group = group)) +
scale_color_manual(name = "instruction", values = instruction) +
scale_fill_manual(name = "instruction", values = instruction) +
labs(x = "\nwave", y = "predicted pattern expression value\n") +
dc_bwsummary(mod_pev_wave)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: dotSTD ~ wave * instruction + (1 + wave * instruction | subjectID)
## Data: filter(data, algorithm == "craving_regulation")
##
## REML criterion at convergence: 46494.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.1474 -0.6296 0.0024 0.6300 4.4455
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subjectID (Intercept) 0.07721 0.2779
## wave2 0.05818 0.2412 -0.37
## instructionregulate 0.10015 0.3165 -0.04 0.03
## wave2:instructionregulate 0.11528 0.3395 0.12 -0.38 -0.38
## Residual 0.72303 0.8503
## Number of obs: 18040, groups: subjectID, 255
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.27366 0.02138 252.53513 -12.801
## wave2 0.08587 0.02445 229.88478 3.513
## instructionregulate 0.89880 0.02645 247.41990 33.983
## wave2:instructionregulate -0.15568 0.03440 221.59634 -4.525
## Pr(>|t|)
## (Intercept) < 0.0000000000000002 ***
## wave2 0.000534 ***
## instructionregulate < 0.0000000000000002 ***
## wave2:instructionregulate 0.00000985 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) wave2 instrc
## wave2 -0.480
## instrctnrgl -0.290 0.245
## wv2:nstrctn 0.264 -0.554 -0.512
Cravings decreased from wave 1 to 2; the effect was stronger when people were looking compared to when they were regulating.
The more the brain looks like a person is regulating, the lower their cravings; this relationships is somewhat stronger (i.e., more negative) when people are looking compared to regulating.
mod_craving = lmer(rating ~ wave * instruction * dotSTD + (1 + wave * instruction | subjectID),
data = filter(data, algorithm == "craving_regulation"))
ggeffects::ggpredict(mod_craving, terms = c("dotSTD", "wave", "instruction")) %>%
data.frame() %>%
ggplot(aes(x, predicted, color = group, fill = group)) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .5, color = NA) +
geom_line(aes(group = group)) +
facet_grid(~facet) +
scale_color_manual(name = "wave", values = algorithm) +
scale_fill_manual(name = "wave", values = algorithm) +
labs(x = "\npredicted pattern expression value", y = "predicted craving rating\n") +
dc_bwggeffects::ggpredict(mod_craving, terms = c("wave", "dotSTD[-1,0,1]", "instruction")) %>%
data.frame() %>%
ggplot(aes(x, predicted, color = group, fill = group)) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .5, color = NA) +
geom_line(aes(group = group)) +
facet_grid(~facet) +
scale_color_manual(name = "pattern expression", values = algorithm) +
scale_fill_manual(name = "pattern expression", values = algorithm) +
labs(x = "\nwave", y = "predicted craving rating\n") +
dc_bwsummary(mod_craving)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rating ~ wave * instruction * dotSTD + (1 + wave * instruction |
## subjectID)
## Data: filter(data, algorithm == "craving_regulation")
##
## REML criterion at convergence: 37063.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8913 -0.6678 0.0148 0.6644 4.0521
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subjectID (Intercept) 0.2142 0.4628
## wave2 0.3191 0.5649 -0.17
## instructionregulate 0.2775 0.5268 -0.45 -0.06
## wave2:instructionregulate 0.1752 0.4186 -0.11 -0.58 -0.28
## Residual 0.4505 0.6712
## Number of obs: 17148, groups: subjectID, 253
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 3.25803 0.03132 253.04648 104.027
## wave2 -0.60967 0.04175 216.05262 -14.604
## instructionregulate -0.92645 0.03723 272.05521 -24.885
## dotSTD -0.03429 0.01213 16012.60486 -2.827
## wave2:instructionregulate 0.10576 0.03732 268.45316 2.834
## wave2:dotSTD -0.01295 0.01764 16177.48242 -0.734
## instructionregulate:dotSTD 0.02846 0.01677 15853.36189 1.697
## wave2:instructionregulate:dotSTD 0.02790 0.02411 14963.52335 1.157
## Pr(>|t|)
## (Intercept) < 0.0000000000000002 ***
## wave2 < 0.0000000000000002 ***
## instructionregulate < 0.0000000000000002 ***
## dotSTD 0.00471 **
## wave2:instructionregulate 0.00494 **
## wave2:dotSTD 0.46287
## instructionregulate:dotSTD 0.08970 .
## wave2:instructionregulate:dotSTD 0.24715
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) wave2 instrc dotSTD wv2:ns w2:STD in:STD
## wave2 -0.235
## instrctnrgl -0.475 0.036
## dotSTD 0.108 -0.082 -0.093
## wv2:nstrctn 0.034 -0.580 -0.382 0.094
## wave2:dtSTD -0.074 0.099 0.064 -0.683 -0.115
## instrct:STD -0.077 0.058 -0.070 -0.718 0.069 0.489
## wv2:nst:STD 0.053 -0.071 0.048 0.493 -0.103 -0.723 -0.685
Pattern expression of the craving signature decreased from wave 1 to 2 (i.e., peoples’ brains looked less like they were craving).
Regulation is associated with lower craving signature expression across waves.
mod_pev_wave = lmer(dotSTD ~ wave * instruction + (1 + wave * instruction | subjectID),
data = filter(data, algorithm == "craving"))
ggeffects::ggpredict(mod_pev_wave, terms = c("wave", "instruction")) %>%
data.frame() %>%
ggplot(aes(x, predicted, color = group, fill = group)) +
geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
geom_line(aes(group = group)) +
scale_color_manual(name = "instruction", values = instruction) +
scale_fill_manual(name = "instruction", values = instruction) +
labs(x = "\nwave", y = "predicted pattern expression value\n") +
dc_bwsummary(mod_pev_wave)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: dotSTD ~ wave * instruction + (1 + wave * instruction | subjectID)
## Data: filter(data, algorithm == "craving")
##
## REML criterion at convergence: 49649.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.2765 -0.6227 0.0061 0.6223 6.7929
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subjectID (Intercept) 0.08926 0.2988
## wave2 0.08692 0.2948 -0.44
## instructionregulate 0.04248 0.2061 0.15 -0.05
## wave2:instructionregulate 0.04861 0.2205 0.22 -0.57 -0.15
## Residual 0.87165 0.9336
## Number of obs: 18040, groups: subjectID, 255
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 0.23606 0.02315 249.22640 10.199
## wave2 -0.05843 0.02808 235.86125 -2.081
## instructionregulate -0.13192 0.02299 249.39700 -5.737
## wave2:instructionregulate -0.01078 0.03185 236.89089 -0.338
## Pr(>|t|)
## (Intercept) < 0.0000000000000002 ***
## wave2 0.0386 *
## instructionregulate 0.0000000278 ***
## wave2:instructionregulate 0.7353
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) wave2 instrc
## wave2 -0.520
## instrctnrgl -0.271 0.259
## wv2:nstrctn 0.322 -0.620 -0.531
The more the brain looks like a person is craving, the higher their cravings (though this relationship was weak). This relationships is somewhat stronger when people are regulating compared to looking.
The relationship between craving signature pattern expression and craving was somewhat stronger at wave 2 compared to wave 1.
mod_craving = lmer(rating ~ wave * instruction * dotSTD + (1 + wave * instruction | subjectID),
data = filter(data, algorithm == "craving"))
ggeffects::ggpredict(mod_craving, terms = c("dotSTD", "wave", "instruction")) %>%
data.frame() %>%
ggplot(aes(x, predicted, color = group, fill = group)) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .5, color = NA) +
geom_line(aes(group = group)) +
facet_grid(~facet) +
scale_color_manual(name = "wave", values = algorithm) +
scale_fill_manual(name = "wave", values = algorithm) +
labs(x = "\npredicted pattern expression value", y = "predicted craving rating\n") +
dc_bwsummary(mod_craving)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rating ~ wave * instruction * dotSTD + (1 + wave * instruction |
## subjectID)
## Data: filter(data, algorithm == "craving")
##
## REML criterion at convergence: 37023.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9327 -0.6625 0.0156 0.6661 3.9978
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subjectID (Intercept) 0.2147 0.4634
## wave2 0.3244 0.5696 -0.17
## instructionregulate 0.2793 0.5285 -0.45 -0.05
## wave2:instructionregulate 0.1797 0.4239 -0.09 -0.59 -0.29
## Residual 0.4492 0.6702
## Number of obs: 17148, groups: subjectID, 253
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 3.26349 0.03128 250.60545 104.342
## wave2 -0.61426 0.04196 215.24030 -14.639
## instructionregulate -0.94102 0.03655 250.48483 -25.743
## dotSTD 0.01721 0.01127 16055.19923 1.527
## wave2:instructionregulate 0.12236 0.03612 229.32164 3.388
## wave2:dotSTD 0.02940 0.01586 16328.59675 1.854
## instructionregulate:dotSTD 0.03003 0.01578 15837.07175 1.903
## wave2:instructionregulate:dotSTD -0.02092 0.02197 15395.51824 -0.952
## Pr(>|t|)
## (Intercept) < 0.0000000000000002 ***
## wave2 < 0.0000000000000002 ***
## instructionregulate < 0.0000000000000002 ***
## dotSTD 0.126901
## wave2:instructionregulate 0.000829 ***
## wave2:dotSTD 0.063818 .
## instructionregulate:dotSTD 0.057023 .
## wave2:instructionregulate:dotSTD 0.341120
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) wave2 instrc dotSTD wv2:ns w2:STD in:STD
## wave2 -0.235
## instrctnrgl -0.483 0.040
## dotSTD -0.083 0.064 0.070
## wv2:nstrctn 0.040 -0.608 -0.368 -0.074
## wave2:dtSTD 0.058 -0.080 -0.049 -0.708 0.092
## instrct:STD 0.058 -0.045 -0.074 -0.708 0.076 0.499
## wv2:nst:STD -0.041 0.056 0.051 0.503 -0.087 -0.713 -0.709
# source functions
library(boot)
source("indirectMLM.R")
# create craving signature dataframe
data_med_craving = data %>%
filter(algorithm == "craving") %>%
mutate(wave = ifelse(wave == "1", 0,
ifelse(wave == "2", 1, NA)),
dotSTD = dotSTD[,]) %>%
select(subjectID, wave, trial, instruction, rating, dotSTD) %>%
data.frame() %>%
na.omit()
# create craving regulation signature dataframe
data_med_regulation = data %>%
filter(algorithm == "craving_regulation") %>%
mutate(wave = ifelse(wave == "1", 0,
ifelse(wave == "2", 1, NA)),
dotSTD = dotSTD[,]) %>%
select(subjectID, wave, trial, instruction, rating, dotSTD) %>%
data.frame() %>%
na.omit()
# define variables
x_var = "wave"
y_var = "rating"
m_var = "dotSTD"model_name = "mediation_craving"
data_med = data_med_craving
if (file.exists(sprintf("models/model_%s.RDS", model_name))) {
assign(get("model_name"), readRDS(sprintf("models/model_%s.RDS", model_name)))
} else {
assign(get("model_name"), boot(data = data_med, statistic = indirect.mlm, R = 500,
y = y_var, x = x_var, mediator = m_var, group.id = "subjectID",
between.m = F, uncentered.x = F))
saveRDS(eval(parse(text = model_name)), sprintf("models/model_%s.RDS", model_name))
}
indirect.mlm.summary(get(model_name))## #### Population Covariance ####
## Covariance of Random Slopes a and b: 0 [-0.001, 0.004]
##
##
## #### Indirect Effects ####
## # Within-subject Effects
## Unbiased Estimate of Within-subjects Indirect Effect: -0.005 [-0.007, 0]
## Biased Estimate of Within-subjects Indirect Effect: -0.005 [-0.008, -0.002]
## Bias in Within-subjects Indirect Effect: 0 [0, 0.004]
##
##
## #### Total Effect ####
## Unbiased Estimate of Total Effect: -0.545 [-0.571, -0.517]
## Biased Total Effect of X on Y (c path): -0.545 [-0.571, -0.516]
## Bias in Total Effect: 0 [0, 0.003]
##
##
## #### Direct Effects ####
## Direct Effect of Predictor on Dependent Variable (c' path): -0.541 [-0.567, -0.513]
## Within-subjects Effect of Predictor on Mediator (a path for group-mean centered predictor): -0.064 [-0.09, -0.033]
## Within-subjects Effect of Mediator on Dependent Variable (b path for group-mean centered mediator): 0.082 [0.067, 0.095]
model_name = "mediation_craving_look"
data_med = filter(data_med_craving, instruction == "look")
if (file.exists(sprintf("models/model_%s.RDS", model_name))) {
assign(get("model_name"), readRDS(sprintf("models/model_%s.RDS", model_name)))
} else {
assign(get("model_name"), boot(data = data_med, statistic = indirect.mlm, R = 500,
y = y_var, x = x_var, mediator = m_var, group.id = "subjectID",
between.m = F, uncentered.x = F))
saveRDS(eval(parse(text = model_name)), sprintf("models/model_%s.RDS", model_name))
}
indirect.mlm.summary(get(model_name))## #### Population Covariance ####
## Covariance of Random Slopes a and b: 0 [-0.004, 0.003]
##
##
## #### Indirect Effects ####
## # Within-subject Effects
## Unbiased Estimate of Within-subjects Indirect Effect: -0.002 [-0.006, 0.001]
## Biased Estimate of Within-subjects Indirect Effect: -0.002 [-0.004, 0]
## Bias in Within-subjects Indirect Effect: 0 [0, 0.004]
##
##
## #### Total Effect ####
## Unbiased Estimate of Total Effect: -0.609 [-0.637, -0.576]
## Biased Total Effect of X on Y (c path): -0.611 [-0.639, -0.579]
## Bias in Total Effect: 0.002 [0, 0.005]
##
##
## #### Direct Effects ####
## Direct Effect of Predictor on Dependent Variable (c' path): -0.607 [-0.636, -0.574]
## Within-subjects Effect of Predictor on Mediator (a path for group-mean centered predictor): -0.053 [-0.092, -0.014]
## Within-subjects Effect of Mediator on Dependent Variable (b path for group-mean centered mediator): 0.032 [0.014, 0.05]
model_name = "mediation_craving_regulate"
data_med = filter(data_med_craving, instruction == "regulate")
if (file.exists(sprintf("models/model_%s.RDS", model_name))) {
assign(get("model_name"), readRDS(sprintf("models/model_%s.RDS", model_name)))
} else {
assign(get("model_name"), boot(data = data_med, statistic = indirect.mlm, R = 500,
y = y_var, x = x_var, mediator = m_var, group.id = "subjectID",
between.m = F, uncentered.x = F))
saveRDS(eval(parse(text = model_name)), sprintf("models/model_%s.RDS", model_name))
}
indirect.mlm.summary(get(model_name))## #### Population Covariance ####
## Covariance of Random Slopes a and b: 0 [-0.004, 0.001]
##
##
## #### Indirect Effects ####
## # Within-subject Effects
## Unbiased Estimate of Within-subjects Indirect Effect: -0.004 [-0.009, -0.002]
## Biased Estimate of Within-subjects Indirect Effect: -0.004 [-0.006, -0.002]
## Bias in Within-subjects Indirect Effect: 0 [0, 0.004]
##
##
## #### Total Effect ####
## Unbiased Estimate of Total Effect: -0.494 [-0.522, -0.468]
## Biased Total Effect of X on Y (c path): -0.495 [-0.523, -0.469]
## Bias in Total Effect: 0.001 [0, 0.004]
##
##
## #### Direct Effects ####
## Direct Effect of Predictor on Dependent Variable (c' path): -0.49 [-0.518, -0.462]
## Within-subjects Effect of Predictor on Mediator (a path for group-mean centered predictor): -0.073 [-0.118, -0.034]
## Within-subjects Effect of Mediator on Dependent Variable (b path for group-mean centered mediator): 0.053 [0.038, 0.067]
model_name = "mediation_regulation"
data_med = data_med_regulation
if (file.exists(sprintf("models/model_%s.RDS", model_name))) {
assign(get("model_name"), readRDS(sprintf("models/model_%s.RDS", model_name)))
} else {
assign(get("model_name"), boot(data = data_med, statistic = indirect.mlm, R = 500,
y = y_var, x = x_var, mediator = m_var, group.id = "subjectID",
between.m = F, uncentered.x = F))
saveRDS(eval(parse(text = model_name)), sprintf("models/model_%s.RDS", model_name))
}
indirect.mlm.summary(get(model_name))## #### Population Covariance ####
## Covariance of Random Slopes a and b: 0 [-0.005, 0.002]
##
##
## #### Indirect Effects ####
## # Within-subject Effects
## Unbiased Estimate of Within-subjects Indirect Effect: -0.002 [-0.008, 0.005]
## Biased Estimate of Within-subjects Indirect Effect: -0.002 [-0.005, 0.005]
## Bias in Within-subjects Indirect Effect: 0 [0, 0.005]
##
##
## #### Total Effect ####
## Unbiased Estimate of Total Effect: -0.551 [-0.568, -0.531]
## Biased Total Effect of X on Y (c path): -0.545 [-0.561, -0.524]
## Bias in Total Effect: 0.006 [0.002, 0.007]
##
##
## #### Direct Effects ####
## Direct Effect of Predictor on Dependent Variable (c' path): -0.549 [-0.564, -0.533]
## Within-subjects Effect of Predictor on Mediator (a path for group-mean centered predictor): 0.008 [-0.02, 0.022]
## Within-subjects Effect of Mediator on Dependent Variable (b path for group-mean centered mediator): -0.234 [-0.243, -0.227]
model_name = "mediation_regulation_look"
data_med = filter(data_med_regulation, instruction == "look")
if (file.exists(sprintf("models/model_%s.RDS", model_name))) {
assign(get("model_name"), readRDS(sprintf("models/model_%s.RDS", model_name)))
} else {
assign(get("model_name"), boot(data = data_med, statistic = indirect.mlm, R = 500,
y = y_var, x = x_var, mediator = m_var, group.id = "subjectID",
between.m = F, uncentered.x = F))
saveRDS(eval(parse(text = model_name)), sprintf("models/model_%s.RDS", model_name))
}
indirect.mlm.summary(get(model_name))## #### Population Covariance ####
## Covariance of Random Slopes a and b: -0.001 [-0.005, 0.002]
##
##
## #### Indirect Effects ####
## # Within-subject Effects
## Unbiased Estimate of Within-subjects Indirect Effect: -0.004 [-0.009, -0.001]
## Biased Estimate of Within-subjects Indirect Effect: -0.003 [-0.006, -0.002]
## Bias in Within-subjects Indirect Effect: 0.001 [0, 0.005]
##
##
## #### Total Effect ####
## Unbiased Estimate of Total Effect: -0.61 [-0.637, -0.572]
## Biased Total Effect of X on Y (c path): -0.611 [-0.637, -0.573]
## Bias in Total Effect: 0.001 [0, 0.004]
##
##
## #### Direct Effects ####
## Direct Effect of Predictor on Dependent Variable (c' path): -0.606 [-0.634, -0.567]
## Within-subjects Effect of Predictor on Mediator (a path for group-mean centered predictor): 0.086 [0.053, 0.124]
## Within-subjects Effect of Mediator on Dependent Variable (b path for group-mean centered mediator): -0.041 [-0.061, -0.022]
model_name = "mediation_regulation_regulate"
data_med = filter(data_med_regulation, instruction == "regulate")
if (file.exists(sprintf("models/model_%s.RDS", model_name))) {
assign(get("model_name"), readRDS(sprintf("models/model_%s.RDS", model_name)))
} else {
assign(get("model_name"), boot(data = data_med, statistic = indirect.mlm, R = 500,
y = y_var, x = x_var, mediator = m_var, group.id = "subjectID",
between.m = F, uncentered.x = F))
saveRDS(eval(parse(text = model_name)), sprintf("models/model_%s.RDS", model_name))
}
indirect.mlm.summary(get(model_name))## #### Population Covariance ####
## Covariance of Random Slopes a and b: 0 [-0.004, 0.003]
##
##
## #### Indirect Effects ####
## # Within-subject Effects
## Unbiased Estimate of Within-subjects Indirect Effect: 0 [-0.004, 0.003]
## Biased Estimate of Within-subjects Indirect Effect: 0 [-0.001, 0.001]
## Bias in Within-subjects Indirect Effect: 0 [0, 0.004]
##
##
## #### Total Effect ####
## Unbiased Estimate of Total Effect: -0.495 [-0.521, -0.466]
## Biased Total Effect of X on Y (c path): -0.495 [-0.522, -0.467]
## Bias in Total Effect: 0 [0, 0.004]
##
##
## #### Direct Effects ####
## Direct Effect of Predictor on Dependent Variable (c' path): -0.494 [-0.52, -0.467]
## Within-subjects Effect of Predictor on Mediator (a path for group-mean centered predictor): -0.062 [-0.097, -0.023]
## Within-subjects Effect of Mediator on Dependent Variable (b path for group-mean centered mediator): 0.003 [-0.012, 0.021]